#include "math.h"
#include "random"
#include "array"
using namespace std;
#include "main.h"
#include "CStemCell.h"
#include "System.h"
extern struct IMD _MD;
extern default_random_engine eng;
extern uniform_real_distribution<double> uniform_dist;
extern normal_distribution<double> normal_dist;
extern uniform_int_distribution<int> int_dist;

extern int global_id ;


void CStemCell::init(PARStemCell & init_par)
{
    this->_par = init_par;
    _cellid = global_id ++ ;
    _X[0] = uniform_dist(eng);       // Initial state of the epigenetic state x1
    _X[1] = uniform_dist(eng)*0.1;    // Initial state of the epigenetic state x2
    _ProfQ = 0;               // We assume that all cells at resting phase at the initial state.
    _age = 0;
}


CStemCell::CStemCell()
        : _X(), _cellid(0), _par(PARStemCell()), _ProfQ(false), _age(0.0)
{
    // 默认构造函数实现
    // 初始化成员变量，如果需要的话
}

CStemCell::CStemCell(const CStemCell& other)
        : _X(other._X),  _cellid(other._cellid), _par(other._par), _ProfQ(other._ProfQ), _age(other._age),
          _dynamic(other._dynamic)
{
    // 拷贝构造函数实现
    // std::array和PARStemCell结构体的成员会自动被复制
}

CStemCell::~CStemCell()
{
}

CStemCell& CStemCell::operator=(const CStemCell& other)
{
    if (this != &other) {
        // 移动成员变量
        _X = other._X;
        _dynamic =other._dynamic;
        _cellid = other._cellid;
        _par = other._par;
        _ProfQ = other._ProfQ;
        _age = other._age;
    }
    return *this;
}
CStemCell CStemCell::getNextStemCell(double mutationrate,double varsigma)
{
    CStemCell cell;
    cell._X = this->getNextEpi();
    cell._par = this->getNextPar(mutationrate,varsigma);
    cell._cellid = global_id ++;
    cell._age = 0;
    cell._ProfQ = 0;
    return cell;
}
void CStemCell::printInfo(std::string msg) const
{
    std::cout << "----" << msg << "-----" << std::endl;
    std::cout << "Stem Cell Information:" << std::endl;
    std::cout << " (1) X: ";
    for (size_t i = 0; i < _X.size(); ++i) {
        std::cout << _X[i] << " ";
    }
    std::cout << std::endl;

    std::cout << " (2) Dynamic: ";
    for (size_t i = 0; i < _dynamic.size(); ++i) {
        std::cout << _dynamic[i] << " ";
    }
    std::cout << std::endl;

    std::cout << " (3) Cell ID: " << _cellid << std::endl;

    // 假设 PARStemCell 类有 printInfo 方法
    _par.printAllElements("(4) 基因参数");

    std::cout << " (5) ProfQ: " << (_ProfQ ? "True" : "False") << std::endl;
    std::cout << " (6) Age: " << _age << std::endl;
}
PARStemCell CStemCell::getNextPar(double mutationrate,double varsigma)
{
    PARStemCell newpar;
    newpar = this->_par;
    int pos = int_dist(eng);
    double varofpar;
    if( uniform_dist(eng) < mutationrate)  // Perform gene mutation with change in parameters.
    {
        varofpar = newpar.getPar(pos);
        varofpar = varofpar * exp(varsigma * normal_dist(eng)) / exp((varsigma * varsigma) / 2);
        newpar.getPar(pos) = varofpar;
        newpar.IsMutated = true;
    }
    return(newpar);
}

std::array<double,NUMVAR> CStemCell::getNextEpi()
{
    std::array<double,NUMVAR> array_re;

    for(int i = 0; i < _X.size();i++)
    {
        array_re[i] = this->getOneEpi(i,this->_X[0],this->_X[1]);
    }
    return(array_re);
}

double CStemCell::getOneEpi(int i,  double x1, double x2)
{
    double phi;
    double a, b;
    double z;
    switch(i)
    {
        case 0:
            phi = 0.08 + _par.d1 * pow(_par.alpha1 * x1, 1.8)/(1+pow(_par.alpha1 * x1, 1.8));
            break;
        case 1:
            phi = 0.08 + _par.d2 * pow(_par.alpha2 * x2, 1.8)/(1+pow(_par.alpha2 * x2, 1.8));
            break;
        default:
            phi = 0;
    }
    a = _par.eta  * phi;
    b = _par.eta * (1-phi);
    // 定义beta分布
    double x,y;

    auto g_dist_1 = std::gamma_distribution<double>{a,1};
    x = g_dist_1(eng);
    auto g_dist_2 = std::gamma_distribution<double>{b,1};
    y = g_dist_2(eng);
    z = x/(x+y);
    return(z);
}


void CStemCell::calcDynamic(double u,int N0, PARStemCell & sys_par,double dt)
{
    double beta0;
    double theta;
    double x1 = _X[0],x2 = _X[1];
    beta0 = _par.beta0;
    theta = _par.theta0 * (1.0 + _par.lambda * pow(_par.a4 * x2, 6.0)/(1.0 + pow(_par.a4 * x2, 6.0)));

    _dynamic[0] = dt*beta0 * (1.0/(1.0 + N0/theta)) * ((_par.a1 * x1 + pow(_par.a2 * x1,6.0))/(1+pow(_par.a3 * x1,6.0)));

    double kappa1,kappa2,dist;

    kappa1 = _par.kappa0 * 1.0/(1.0 + pow(_par.b1 * x1, 6.0));
    dist = _par.calcDistance(sys_par);//
    kappa2 = _par.phi0 *(pow(dist,4))/(pow(dist,4)+pow(_par.z0,4));
    _dynamic[1] = dt*(kappa1 + kappa2);
    double g;
    g = 3.0 * pow(x2, u) * pow(1-x2, 1 - u);
    _dynamic[2] =dt*_par.mu0/(1 + _par.c1 * exp(g));
}